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Muller's ratchet is a paradigmatic model for the accumulation of deleterious mutations in a pop- 
ulation of finite size. A click of the ratchet occurs when all individuals with the least number of 
deleterious mutations are lost irreversibly due to a stochastic fluctuation. In spite of the simplicity 
of the model, a quantitative understanding of the process remains an open challenge. In contrast 
to previous works, we here study a Moran model of the ratchet with overlapping generations. 
Employing an approximation which describes the fittest individuals as one class and the rest as 
a second class, we obtain closed analytical expressions of the ratchet rate in the rare clicking 
regime. As a click in this regime is caused by a rare large fluctuation from a metastable state, we 
do not resort to a diffusion approximation but apply an approximation scheme which is especially 
well suited to describe extinction events from metastable states. This method also allows for a 
derivation of expressions for the quasi-stationary distribution of the fittest class. Additionally, we 
confirm numerically that the formulation with overlapping generations leads to the same results 
as the diffusion approximation and the corresponding Wright-Fisher model with non-overlapping 
generations. 



I. INTRODUCTION 



In an asexual population of finite size, weakly deleterious mutations can fix by chance. This phenomenon is due to 
stochastic fluctuations originating from the finiteness of the population, which can lead to a loss of the fittest class of 
individuals. If the genome is assumed to be so long that back mutations are unlikely and can be ignored, the fittest 
class is lost forever and the number of fixed deleterious mutations increases irreversibly. This process has been termed 
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Furthermore it has 



seen thought to a ccount for t he degeneration of non-re combining pa rts of sexually reproducing 
organisms such as the Y-chromosome (Rice 19871 and mitocho ndrial DNA ( Lynch] 1996). Muller's ratchet can also 
be used to explain the absence of long-lived asexual lineages (Lynch et al. 1993). Since in the absence of back 



mutations, mutation-free genomes can only be recreated by recombination between mutation-loaded classes, Muller's 
ratchet provides an appealing explanation for the evolution of sex (Barton and Charlesworth 1998 De Visser and 
Elenal 120071). 



Each time the least-loaded class, i.e. the class with the fewest number of deleterious mutations, is lost, it is said 
that Muller's ratchet has clicked. Since the rate of the ratchet determines the speed of degeneration of the population, 
this quantity is of central interest. In its simplest form the rate of Muller's ratchet depends only on the selection 
coefficient S, the mutation rate U and the size N of the population, where it is assumed that each mutation has the 
same effect so individuals with k mutations have fitness {1~ S)''. In this case the fitness space is equivalent to an axis 
counting the number of deleterious mutations and the population can be organized into discrete classes labeled by the 
number of mutations they carry. The deleterious mutations have the effect of shifting the population to higher values 
of k. Since the fitness of the respective classes is given by (1 — S)'', selection works into the opposite direction. In 
the limit of an infinitely large population these two opposing forces lead to a steady state distribution whose precise 
form was found by Haigh (Haigh 1978). 



If finite populations are considered, however, the calculation of the rate of Muller's ratchet turns out to be an 
intricate problem, despite its simple formulation. The difficulty arises due to the complex interaction of the fluctuation 
of the least-loaded class with the rest of the distribution. A detailed quantitative understanding of the behavior of 
the occupation of the class with the fewest mutations, however, is necessary to determine the mean time to extinction 



of this class, i.e. the inverse of the ratchet rate. Despite of considerable efforts and recent adva nces (^Etheridge et al. 
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ratchet rate remains a challenging open prob. 

In its standard form Muller's ratchet was first quantitatively described by Haigh who analyzed a classical Wright- 
Fisher model of an asexually reproducing population of fixed size N. He pointed out that the most important quantity 
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of the ratchet is the average number of individuals in the least loaded class, no = N ex-p{—U / S) , because fluctuations 
of no ultimately lead to a click of the ratchet. Later it was shown by Jain that the ratchet rate cannot depend only 
on no but rather has to depend on uqS (Jain 2008). If UqS is small, the ratchet clicks frequently and the populations 
behaves like a wave in /c— space propagating towards higher values of k. The traveling wave approach to Muller's 
ratchet was discussed in ( jRouzine et a/., ,2008, ) and provides an appealing quantitative theory for frequently clicking 
ratchets. 

While this regime of Muller's ratchet is relatively well understood, a quantitative understanding of the opposite 
case no 5' ^ 1 of a rarely clicking ratchet is still lacking and has recently attracted a lot of atten tion (Jain, |2008{ 



Neher and Shraiman 2012). In this regime the rate of the ratchet is exponentially small in noS ( Jain[^ 2008 ) and 



extinction of the fittest class occurs as the result of a rare large fluctuation. In contrast to the fast clicking regime 
the distribution of the population equilibrates to a metastable state after each click. A wide-spread approach in this 
regime is therefore to consider only the fittest class and apply a phenomenological model for all the classes nfc, fc > 0, 
with more mutations than the least-loaded class. Such an approach leads to a one-dimensional approximation where 
just the fittest class is taken into account. Generally the rate of the ratchet can then be calculated by means of a 
diffusion approximation as the result of a one-dimensional mean-first passage problem. Recently it was shown how 
this approach can be improved by accounting for the interaction of the fiuctuations of the fittest class with the tail of 
smaller fitness which can lead to a delayed feedback (Neher and Shraiman 2012). 

Up to now a quantitative treatment of Muller's ratchet relied either on Haigh's model or on the corresponding 
diffusion approximation. To our knowledge a Moran formulation with overlapping generations has not been employed 
so far. This is not surprising as in the diffusive limit any quantity should become independent of the respective micro- 
scopic formulation and a Moran formulation of Muller's ratchet is expected to be computationally disadvantageous. 
A Moran formulation, however, can also lead to interesting new approaches to tackle the problem of the ratchet rate 
analytically. 

In the present work we investigate a Moran formulation of Muller's ratchet and show how this model can be 
approximated by a one-dimensional Moran-process in the regime no5 ^ 1 where the ratchet clicks infrequently. We 
show that this model allows for an analytical solution for the ratchet rate which agrees almost perfectly with values 
obtained by numerical simulations of the full ratchet. Furthermore, by employing a recently developed method to 
treat rare large fiuctuations in stochastic population dynamics, we find analytical expressions for the ratchet rate and 
the quasi-stationary distribution of the fittest class in the parameter range U/S < 2 which also agree very well with 
the corresponding results of the full ratchet. Finally, we confirm numerically that the formulation with overlapping 
generations leads to the same results as the diffusion approximation and the corresponding Wright-Fisher model with 
non-overlapping generations. 



II. MODELS AND METHODS 



In the standard formulation of Muller's ratchet, as considered by Haigh (Haigh', '1978'), mutations in a population of 
fixed size N occur at rate U and individuals are classified into different groups according to the number of deleterious 
mutations they carry, k. Each mutation reduces the fitness of the genotype k by an amount S such that the growth 
rate of an individual with k mutations is proportional to (1 — S)''. The ratio of the mutation rate U to the mutation 
effect S, which is denoted by A = U/S, plays a central in the analysis of the ratchet. 

The reproduction model usually employed in the analysis of Muller's ratchet is Wright-Fisher sampling. It consists 
of, at each time step, replacing the whole generation of individuals by a multinomial resampling of the current 
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and the mean fitness W is given by = X]j!Lo(l ~ SyPi- 

Wright-Fisher sampling has the advantage of being very efficient for numerical simulations. The downside of the 
model is, however, that it does not easily allow for analytical methods to be used. Therefore, the corresponding 
diffusion approximation of the microscopic Wright-Fisher formulation is usually used to predict the click rates of the 
ratchet. 

The second widely applied reproduction model in population genetics is the Moran process which in contrast to the 
Wright-Fisher formulation assumes overlapping generations. The Moran process on which we focus in this article is 



amenable to a wider range of analytical methods (at the cost of being slower in numerical simulations) ( Park et al. 
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2010). It is a stochastic process in which at each time step one individual is chosen for reproduction and one for 



removal from the population. The choice of the individual that reproduces is random, but (similarly to the Wright- 
Fisher formulation) weighted by the fitness of the class the individual is chosen from. The probability of removal (or 
death) of an individual is independent of the fitness. Applied to Muller's ratchet this therefore embodies the following 
procedure: An individual with k mutations is chosen according to the abundance and selection preference of the class 
k with weight (1 — S)''nk / '^'^q{1 — Syuj. This individual spawns one offspring with k mutations that can then 
mutate to k + I mutations with probability /II. The probability to mutate is thus 1 — e^^ , which is the same 

as in the Wright-Fisher model. Also, one individual with k mutations is chosen for removal with probability Uk/N 
(this may be the one that reproduced). Since on average every individual is chosen for removal once every N time 
steps, it is natural to define one generation in the Moran model as N time steps. In all figures, the ratchet click times 
are thus expressed in generations. 

Although different on the microscopic scale, both Wright-Fisher and Moran models usually converge to the same 
mcsoscopic diffusion regime when A'^ is large and fitness advantages and mutation rates are of order N~^. In this 
limit, the equation describing the evolution of the population is given by 



dt 



Hk = S{k - k)nk -Uuk + Uuk-i + y/n^ijk 



(2) 



where k = J^k ^'^fc (Neher and Shraiman 2012). The uncorrelated Gaussian white noise ij with {'rik{t)rji{t')) — 
Ski^{t — t') models the stochastic fiuctuations due to the finiteness of the population (genetic drift). In the infinite 
population limit, this equation becomes deterministic and has the steady-state solution fik = Ne~^ A''/k\ (Haigh 
1978). Also, a time dependent solution of the deterministic model has been obtained (Etheridge et al. 2009). In this 



Also, a time dependent solution of the deterministic model has been obtained ( Etheridge et al. 
paper, we solve Eq. (|2| numerically using stochastic Runge-Kutta methods (RoCler 2010). 



III. APPROXIMATIVE ONE-DIMENSIONAL MORAN MODEL OF MULLER'S RATCHET 



A mathematical analysis of the Moran model for Muller's ratchet is complicated and even the formulation of the 



corresponding Markov chain (Ewens 2004) is involved and barely leads to new insights. The important advantage of 
Moran models, however, is that they can be formulated in terms of a master equation (Gardiner 2009). Therefore 



Moran models are analytically tractable even beyond the diffusion approximation, if only two species are considered. 
Thus an appealing approach to the analysis of the rate of Muller's ratchet is to approximate the full ratchet by a 
model consisting only of two species. Since we are interested in the loss of the fittest class with zero mutations a 
natural choice is to consider individuals with zero mutations as one species, and to combine all others in a class which 
contains all individuals with one or more mutations. The constraint of a fixed population size N then leads to a 
one-dimensional model. An illustration of the distribution of individuals in the space of mutations and the reduction 
to a one-dimensional model is given in Fig. [T] 

The master equation for the probability P{n, t) to find n = uq individuals in the fittest class is 
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where u is the mutation rate away from the fittest class and s is the fitness disadvantage when carrying a mutation 
( [Blythe and McKane[ [2007] ) . (H ere and in the following lowercase letters denote the parameter values of the approx- 
imative two-state model, while capital letters denote the parameters of the full ratchet.) The idea of representing all 



classes but the fittest as one class was first introduced in (Waxman and Loewe 2010) for a Wright-Fisher model of 
Muller's ratchet. 

A crucial step in the reduction of the full model of Muller's ratchet to the one-dimensional formulation is the relation 
of the two mutation rates and fitness disadvantages in the respective models. This mapping is not unique and two 
reasonable assumptions have to be invoked to relate the two parameters pairs. A plausible approach is to compare 
the steady state distributions in the infinite population limit of the respective models. For the full ratchet the well- 
known steady state distribution for the probability of an individual to have k mutations is Ps{k) = e~^^^ {U / S)^ /k\. 
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FIG. 1 Illustration of MuUer's ratchet in the space of deleterious mutations. Individuals are grouped into different classes 
depending on the number of mutations k they carry. Mutation (blue arrows) drives the population to higher values of k, while 
selection (yellow arrows) opposes this motion, leading to a quasi-stationary distribution. The two-state approximation amounts 
to putting all mutated individuals in one mutated class (light green box). Both mutation into this class and selection pressure 
operating on it (large arrows) have to be calculated from the original mutation rates and selection strengths. Since the total 
population size is conserved, calculating the distribution of the number of individuals in the two classes reduces to the analysis 
of no and thus to a one-dimensional problem. 



A non-zero steady state of the fittest class in the two-state system can only be obtained in the parameter regime 
u/s < 1 and is given by n* = iV(l — u/s). To relate the parameters we now demand that (i) the mean fitness of the 
full population and (ii) the mean fitness of all individuals carrying a mutation is equal in both models. The mean 
fitness of the full population in the steady state of the full ratchet is J2T=QPi^)i^ ~ ^)'' ~ while the mean fitness 
of all individuals in the two state model is N^^{n* + (1 — s)(l — n*)) = 1 — u. Condition (i) accordingly suggests the 
relation 

u = l-e-^ . (6) 

The mean fitness of all individuals carrying mutations is in the full ratchet model given by ~ — 

— 1^ . In the two state model this corresponds to (1 — s) {N — n*)/N = (1 — s)u/s. Employing condition (ii) 
consequently yields the relation 

1 - e-^ 

We can also introduce the parameter X — u/s which is related to A — U/S according to 

A = 1 - e-^ (8) 

Relation ([s]) shows that the restriction A < 1 of the two-state model does not restrict the range of A. 

Before we present the analytical solution for the ratchet rate of this model, let us shortly discuss the validity of 
the approximation used. To correlate the parameters of the full ratchet and the two state model, we have related 
properties of the equilibrium solution of an infinite population in both models. This certainly makes sense as long as 
the typical time tr that it takes for the population to relax to a metastable state after each click is much smaller than 
the mean time t between two successive clicks. This condition is fulfilled in the case of the slowly clicking ratchet, 
which is the regime we focus on in this work. If the ratchet clicks rapidly the population does not equilibrate after a 
click and relating the parameters based on equilibrium distributions is clearly not valid. 

A consequence of our relation (|6| is that the mutation rates out of the fittest class are equal in both models which 
certainly is a reasonable assumption ( [Waxma n and Loewe 20101. Furthermore our second relation ([t]) entails that 



the number of individuals which are not in the fittest class is the same in the equilibrium states of both models, 
i.e. uq = n* . Consequently the same holds true for the number of individuals in the fittest class. Thus, although the 
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parameter mapping is not unique, it is hard to think of any other relation in the slowly clicking regime as this would 
consequently violate the properties specified above. It is important to keep in mind that the parameter mapping is 
only valid in the rare clicking regime and that other mappings might be more appropriate in the fast clicking regime 
(Waxman and Loewe 2010). 



A. Exact solution and Comparison with Full Moran ratchet 



With the reduction to a two-state problem as given in the previous section, we can now exploit the advantages that 
the Moran formulation offers for analytical calculations. The mean click time of the ratchet is given by the mean first 
time of the population with no mutations reaching zero. It is well known that a solution of such a mean first passage 
time problem in a two-state model can be formally written as a product of the transitions rates Q and ^ and is 



given by ( Gardiner , 2009 1 



Af-l 



(9) 



z=y 



where (f>{n,) — YVk=i '^-(^)/^+(^) ^^'^ the initial number of individuals in the lowest mutation class population. 
This expression can be evaluated for moderate but the number of terms grows quickly with N which makes it 
more and more difficult to evaluate t. 

To compare our analytical results to the full ratchet we have performed extensive numerical simulations of the full 
Moran ratchet using the rules detailed in section We organize our results as follows: The parameters U and S 
are grouped according to the conditions specified below, a nd then N is varied. Since the selection penalty S can 

we group parameters with the same rescale d mutation 



2012) 



be interpreted as a timescale (Neher and Shraiman 
rate A = U/S. Similarly, since {N S)~^ can be interpreted as rescaled variance of the stochastic effects (Neher and 
Shraiman 2012), we also group parameters with the same N S, which is then equivalent to keeping uq S = NSe~-'^ 



fixed. The corresponding two-state parameters are rescaled as given by ^ and ([t]). A comparison of the analytical 
results and the simulations is given in Fig. [2j We observe excellent agreement of the analytic result given by Eq. ^ 
for the two-state model with the simulation of the full ratchet in the slow ratchet regime, where the two-state 
approximation is valid. 
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Exact solution 
A=0.1, re|,S=9.05 
A=0.25, n„5=7.79 
A=0.5, n„S=6.07 
A=1.0, re„S=3.68 
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FIG. 2 Comparison of the analytical expression of the mean click time (i.e. the inverse ratchet rate) for the two-state model and 
numerical simulations of the full ratchet in the slow ratchet regime. Different colors correspond to different effective mutation 
rates A = U/S and different fioS = NSe~^. 
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IV. WKB-APPROXIMATION TO THE RATCHET RATE 



The expression ^ for the mean time to extinction is exact. It gets, however, unwieldy and impractical when larger 
population sizes arc considered. Furthermore it does not allow for any analytical statements about the distribution 
of the frequency of the fittest class. Therefore, we want to gain quantitative insight into the ratchet rate and the 
distribution of the frequency of the fittest class by an approximative treatment of Eq. ([s]) . The most widely applied 
approach certainly is the diffusion approximation from which by standard methods the mean time to extinction (MTE) 
T can be calculated analytically ( (Gardiner, |2009) . The resulting expression usually has to be evaluated numerically. 
While the diffusion approximation provides faithful results in the regime where an extinction event is the outcome of 
a typical fluctuation of the process, it in general may fail to describe the MTE correctly when extinction occurs as 
the result of a rare large fluctuation (Doering et al. 2006 Hanggi et al. 1984 Ovaskainen and Meerson 2010). In 



the rare clicking regime the relaxation time to the metastable state is much shorter than the mean time between the 
clicks and the population equilibrates after each click. It is important to note that in such a scenario the click of the 
ratchet is due to a rare large fluctuation away from the metastable state. 

An approach to the treatment of master equations which is especially well suited to account for rare event statistics 
is the WKB- (Wentzel-Kramers-Brillouin) theory. This approximation scheme which is sometimes referred to as 
the eikonal approximation was first developed for a semi-classical treatment of quantum mechanics and has recently 



attracted a lot of attention in the context of stochastic population dynamics (Assaf et al. 2011 Black and McKane 



2011 Dykman et al. 2008 Schwartz et al. 20111. Similar to the diffusion approximation, it replaces the master 
equation of the Moran process by an analytically tractable equation which in addition allows for a mathematically 
controlled approximation in terms of powers of the inverse population size. Recently the WKB-approximation has 



and later considerably extended and generalized in (Assaf and Meerson 2010). 
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The basic idea relies on the fact that the process can be characterized by a metastable state around which the 
frequency of the fittest class resides. After a long average time r the fittest class is eventually lost and MuUer's 
ratchet clicks. For the approach to work two crucial assumptions have to be made. First, the population size has to 
be finite and not too small, i.e. ^ 1. Second the typical relaxation time tr to the metastable state should be much 
shorter than the MTE, i.e. tr <^ r. We note that here this condition has to hold anyway in order for the two state 
approximation to be meaningful. It can be shown that the metastable state, which is sharply peaked around n*, is 
encoded in the first excited eigenvector n(n) of the master equation ^ which has not decayed at a time t > tr (Assaf 



and Meerson 2010). Thus the shape of the PDF of the metastable state, which is referred to as the quasi-stationary 
distribution (QSD), is given by n(n). Furthermore, the decay rate of this distribution, i .e. the ratchet rate r~^, i s 
determined by the first non-zero eigenvalue of the master equation ([3| . As was shown in ( Assaf and Meerson 2006 1 , 
the decay of the QSD for times t ^ tr can therefore be obtained as 



P{n > 0,t) = n(n)e~*/^. 
Accordingly, the click probability distribution behaves as 

P(0,i)==l-e-*/^. 
Using Eqs.(|3|, ([To]) and (11) the click rate is given by 

T-^ = r_ (n = l)n(n = 1) 



(10) 



(11) 



(12) 



which is just to the probability flux into the absorbing state n = 0. 

In following we present an approximative approach to calculate the QSD n(n) based on a WKB-type approximation. 
Inserting (10) into ([s]) one obtains after introducing x = n/N 



^^^^ -t^(x+^\Ti(x+^ \ + t+\x 



(13) 



where n(ri) = Ii{Nx) = 7r(x) and T±{n) = T±{Nx) — t±{x). Since we consider the rare-clicking regime of the ratchet, 
the term on the left-hand side is exponentially small in N and we can neglect it. The resulting quasi-stationary master 
equations reads 



= i_ 
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Now we are ready to employ the WKB approach by expressing the sohition of this equation by the ansatz (Dykman 



et al. 1994) 



Tr{x) = A{x)exp{~N[Sa{x)+0{N-^)]) 
= Cexp{-NSQ(x)~Si{x) + 0{N~^)) 



(15) 



where both Sq{x) and Si{x) are assumed to be of order unity and C is a normalization constant. Inserting this ansatz 
into (14), expanding So{x + l/N) around x to first order and neglecting terms of order 0{1/N), we obtain in leading 



orde 

= t+{x) (e^«(^) " + ^-('^) (^■^"^''■' " 1 
where S'q{x) — -^So{x). The solution of this equation is given by 



So{x) = 



t-iy)\ 



dy. 



(16) 



(17) 



After insertion of So{x) into the ansatz ( 15 ) the lowest order solution for the QSD is obtained up to the normalization 
constant C. To determine C one exploits the fact that the QSD is strongly peaked around x* = n* /N and then 
assumes it to be of Gaussian shape centered at x* which is normalized to unity. Around the maximum a; ~ a;* 
this leads to an approximation of the QSD by 7r(a;) ~ Ce"^'^"'^ )-iN/2)Sg {x ){x-x ) -^^j^ose normalization yields 



C 



-NSo{x') 



Hence in leading order we obtain for the QSD 



tt{x) 



^NlSoix')~So(x)] 



Using this expression of the QSD we can calculate the leading order behavior of the click rate 

^-1 ^ ^-NlSo(x')~Som 



(18) 



(19) 



where we have used that tt{N^^) ^ 7r(0) and t^{N^^) ^ t'_{0)/N for large N. In leading order we thus obtain the 
anticipated exponential dependence of the ratchet on N in the rare clicking regime. These results are valid as long 
as N[So{x*) — So{0) » 1 because the WKB-ansatz requires the ratchet rate to be exponentially small. Furthermore 
the normalization procedure can be expected to fail if the metastable state is close to boundary x* ~ 1 because the 
Gaussian approximation does not hold anymore. 

The next order C'(l)-corrections of the WKB-approximation are obtained by expanding So{x + l/N) to second 
order and Si{x + l/N) to first order around x and provide the pre-factor of the QSD and the ratchet rate. The 
calculation of the sub-leading corrections is more involved and shall not be carried out in detail here. The crucial 
step in the calculation is to note that the WKB-solution in leading order is not valid close to the absorbing state at 
X — 0. Therefore the WKB solution has to be matched with an exact recursion solution of quasi-stationary master 
equation ( 14 ) . A detailed account of this method is given in ( Assaf and Meerson 2010). Following the steps in ( Assaf 
and Meerson, 2010) we obtain for the QSD 



^ ^ y/S"{x*)t+{x*) ^NlSoixn-Solx)] 

^2t: Nt+{x)t^{x) 



The WKB solution for the inverse of the ratchet rate is given by 

^/2'KNR 



^/S^)U{x*){R^l) 



^N[Sa{x*)^Som 



(20) 



(21) 



with i? = i'_(0)A'^(0) = (u- l)/(s- 1). 

Inserting the respective transition rates, we obtain for the mean time to extinction of the fittest class, i.e. the inverse 
of the ratchet rate 



Ns{l - u){N{s - 1) - .s) u/s{l - s)^^T^ {u/s{l ~ u)) 



-AT 



^u/(27riV(l - s))(s - uY{{N - l)(s - 1) - u) 



(22) 
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This expression provides an exact result to order N~^. In Fig.[3|we have compared this result for different parameters 
to the numerical results of the full ratchet. The WKB approximation of the mean time time to extinction in the two 
state model agrees in the range A > 0.5 corresponding A < 1 almost perfectly with the numerical results of the full 
ratchet. While the WKB-prediction is still quite good for A = 2 it starts to deviate for increasing values of A. The 
parameter range in which the WKB-theory works thus is slightly more restricted than in the two-state model. This 
can be explained by noting that for A > 1 the two-state approximation is still valid if the ratchet operates in the rare 
clicking regime, i.e. if NS is chosen to be large enough, see Fig. [2j The WKB-theory on the other hand breaks down 
if X* = e"^ is close to the absorbing state at a; = independent of NS. 
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FIG. 3 Comparison of the WKB-solution (221 for the inverse ratchet rate and the numerical simulations of the full ratchet. 
Different colors correspond to different effective mutation rates K = U/S and different noS = NSe~^ 



To gain a deeper understanding of the WKB-solution one can simplify the unwieldy expression (22 1 for N ^ 
Expanding in N and keeping only the leading order term, we obtain the approximation 



1. 



27rA^ 



(s - uf 



„Ar(s-u(l-log(n/s))) 



(23) 



which is almost indistinguishable from the WKB-solution (22| for > 100. A comparison of the exact solution 
([9|, the WKB-solution (22) and the approximation (23) is provided in the appendix. It is important to note that 



expression (23) exhibits the same scaling behavior in n*s as the one found by Jain in the rare clicking regime (Jain 
2008 ) . To see this we rewrite this expression as 



A^V2^/3(A)s~5(„*)-2e' 



n* sa.{\) 



(24) 



where a(A) = (1 — A[l — lnA])/(l — A) and /3(A) = ^A(l — A'^) . Our solution exhibits A-dependent functions in 
the exponent and the pre-factor which is in contrast to the result of Jain where these factors have to be replaced by 
a constant b. Since a(A) G [0, 1] and /3(A) G [0, 1] the values of both functions are close to the values between b=0.5 
and b=0.6 which were ad hoc chosen for this constant. A dependence on A of the factor b was also observed recently 
by Neher and Shraiman ( |Neher and Shraiman[ 2012). 

The WKB-theory not only yields results for the ratchet rate but is also capable of describing the frequency distri- 
bution of the fittest class in the metastable state, i.e. the QSD, because the parameters of the two-state model were 
chosen such that the size of the fittest classes match in both models. One can therefore expect that also the QSD 
of the fittest class is approximately the same in both models. In Fig. [4] we have compared the numerically obtained 
size of the fittest class in the full ratchet model with the WKB-solution ( p0| ) and observe a striking agreement. As 
anticipated the WKB-theory starts to deviate if the deterministic fixed point x* is close to the absorbing point at 
x=Q and if x* ~ 1. 
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FIG. 4 Comparison of the WKB-solution (201 for the QSD ■k{x) (blue line) with the distribution of the fittest class of the 
full ratchet obtained by numerical simulations for 10^ realizations (green histogram) at 0.1% of the respective click times. At 
this time, the distribution of the class with the lowest number of mutations has already relaxed to the quasi-stationary state, 
while in almost no realization a click has already occurred. The parameters used are (a) A'^=100, A=0.1, (b) A'^=100, A=0.2, 
(c) A'^=200, A=0.5, (d) Af=500, A=l, and in all cases S^O.l. The analytic solution of the two-state model fits the numerical 
distribution obtained for the full ratchet very well. Deviations occur when the fixed point of the deterministic solution, x* , 
begins to approach absorbing point at a:: = 0, which is where the WKB approximation is expected to break down. 



V. COMPARISON OF MORAN, WRIGHT-FISHER AND DIFFUSION MODELS 

Previous studies of Muller's ratchet have mostly considered the diffusion equation ^ to approach MuUer's ratchet 
analytically, while numerical simulations have relied on Haigh's model with non-overlapping generations using Wright- 
Fisher sampling, Eq. ([T]). Since in this work we have used a Moran formulation of the ratchet with overlapping 
generations, we now proceed to compare our results with the Wright-Fisher and diffusion formulations. We note 
that some care has to be taken to ensure that the diffusive limit of the Wright-Fisher model has the same diffusion 



constant as the corresponding Moran formulation, since these usually differ by a factor of two (Ewens 2004). Since 
fluctuations scale with N~^/^, one possibility to take this into account is to consider the Wright-Fisher model with 
N/2 individuals, which is what we do in the simulations presented below. 

We have performed numerical simulations of the Wright-Fisher model and have numerically integrated the diffusion 
equation ^ using stochastic Runge-Kutta methods. To compare the three different approaches, the click times 
averaged over 1000 realizations for each model for different values of A and NS similarly to the previous sections are 
presented in Fig. [5] We observe excellent agreement of the two macroscopic models and the diffusive description for 
slow and fast ratchets. 



VI. CONCLUSIONS AND DISCUSSION 

Muller's ratchet is a model where a single rare large fluctuation influences the fate of the whole population. Therefore 
a quantitative treatment of this model needs to correctly account for such rare fluctuations which renders the analysis 
highly non-trivial despite its simple formulation. In this article we have considered a Moran formulation of Muller's 
ratchet. When the ratchet clicks infrequently the population relaxes to a metastable state after each click event. 
For this regime, where a click can be attributed to a rare large event, we have proposed an approximative model 
which can be described by a one-dimensional master equation. The exact solution of this model for the average time 
between successive clicks agrees almost perfectly with the numerical results of the full ratchet model. To gain further 
insight we furthermore treated the proposed model perturbatively, employing a recently developed scheme which is 
particularly well suited to the description of rare large fluctuations. We obtained a closed-form expression for the 
inverse ratchet rate which is in excellent accordance with the results obtained by numerical simulations of the full 
ratchet. It is worth pointing out that the employed WKB-theory is a non-diffusive approximation in contrast to the 
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FIG. 5 Comparison of the average click times t of MuUer's ratchet for the Moran model (circles), Wright-Fisher model (squares) 
and the diffusion limit of the two models (crosses). Same sets of parameters U and S are shown in the same color. We observe 
perfect coincidence of all three models for both slowly and fast clicking ratchets. 



diffusive approaches applied in most other vi^orks on the rate of Muller's ratchet. Additionally the WKB approach 
allows for an analytical calculation of the frequency distribution of the fittest class in the quasi-stationary state. Given 
that the two approximations, namely the reduction to two classes and the WKB-approach, are applied, the agreement 
of the analytical result with the numerical simulation of the full ratchet model is striking. This finding, in turn, 
strongly supports the validity of the reduced two-state model. To our knowledge this is the first time that analytical 
expressions for the quasi-stationary distribution of the fittest class could be calculated. The expression (20 1 could 
also be of interest for future experiments where the average time between successive ratchet clicks is large because the 
ratchet relaxes very quickly to its metastable state. As we have shown this relaxation time is only a fraction of the 
average inter-click-time and therefore might be experimentally more accessible in the slow ratchet regime. Finally, we 
have confirmed that our results are in excellent agreement with the diffusion approximation and the corresponding 
Wright-Fisher formulation of the ratchet. 



VII. APPENDIX 



A. WKB-approxImation and exact solution of the two-state model 

We compare the exact solution of the click times in the two-state model given by Eq. ([9| to the analytic approxi- 
mation calculated from the WKB-approximation, Eq. (22), and to the approximation of the WKB-result for large N, 
Eq. ([23|, in Fig. [6) The WKB results are in excellent agreement with the exact solution for the rare clicking regime 
when the equilibrium point of the deterministic equation, x* = is sufficiently far away from 0. For large A, the 
WKB approximation begins to deviate from the exact solution. 

For the approximation of the WKB solution for large N, which yields the simple expression Eq. ( 23 1 , we observe 
excellent agreement with the full WKB solution when TV > 100. 
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